##### PRELIMS ##########################################################

rm(list = ls())

library(haven)
library(sp) 
library(sf) 
library(rgdal)
library(raster)
library(readxl) 
library(dplyr)
library(classInt)
library(RColorBrewer) 
library(tmap) 
library(Hmisc) 
library(GISTools) 
library(viridis) 
library(viridisLite) 

tmap_mode("plot") 

# main directory here
dir <- paste0("path_to_main_directory_here/")
dir <- paste0("~/Dropbox/2nd year paper/restat replication files/")

# set directory
setwd(dir)

##### PREPARE DATA ##########################################################

# load shape file
ProvShape <- readOGR(paste0(dir,"orig/additional data/shape"), "Provincias_ETRS89_30N")
# remove Canary Islands and Ceuta and Melilla
ProvShape <- ProvShape[-c(35, 38, 51, 52), ]
# get numerical data
dataformaps_prov <- read_dta(paste0(dir,"data/int/GeoDistSUB.dta"))
# merge
ProvShape <- merge(x = ProvShape, y = dataformaps_prov, by.x = "Codigo",  by.y = "province", all.x = TRUE)

##### GENERATE MAP ##########################################################

map <- tm_shape(ProvShape) +
  tm_fill(col = "ShrPplSUBProv",
          title = "Share in Province (%)",
          breaks = c(0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8),
          legend.show = TRUE) +
  tm_shape(ProvShape) +
  tm_borders("black", lwd = .1) +
  tm_layout(legend.title.size = 0.8, legend.text.size = 0.7,
            legend.bg.color = "white",
            frame = FALSE) +
  tm_legend(position = c("right", "bottom"))

# Save map
tmap_save(tm = map, filename = paste0(dir,"figs/SUBGeoDist.pdf"), width = 7, height = 5)

##### CLOSE ##########################################################
rm(list = ls())